from docplex.mp.model import Model
import time

def append_to_file(filename, content):
    with open(filename, 'a') as file:
        file.seek(0, 2)
        file.write(content)
        file.write('\n')  # 可选择性地添加换行符

def solver(n, verticesNum, params: dict):
    M = 1000
    length = verticesNum + 1
    model = Model(name='LP')

    # x_{ij}^{mn}
    # c
    x = {}
    x_ret = {}
    for i in range(length):
        x[i] = {}
        x_ret[i] = {}
        for j in range(length):
            x[i][j] = {}
            x_ret[i][j] = {}
            for k in range(length):
                x[i][j][k] = {}
                x_ret[i][j][k] = {}
                for l in range(length):
                    x[i][j][k][l] = model.binary_var(
                        name=f'x_{i}_{j}_{k}_{l}')
                    x_ret[i][j][k][l] = 0

    # sd对mn，表示业务量
    sd = []
    for i in range(1, length):
        for j in range(1, length):
            if i != j:
                sd.append((i, j))

    # p_{ij}^{mnvt}
    # 修复路径所占的比例，mn索引到对应的SD对中，vt索引到SD路由方案上的某条边，ij索引到被分流的边
    p = {}
    p_ret = {}
    for i in range(length):
        p[i] = {}
        p_ret[i] = {}
        for j in range(length):
            p[i][j] = {}
            p_ret[i][j] = {}
            for k in range(length):
                p[i][j][k] = {}
                p_ret[i][j][k] = {}
                for l in range(length):
                    p[i][j][k][l] = {}
                    p_ret[i][j][k][l] = {}
                    for m in range(length):
                        p[i][j][k][l][m] = []
                        p_ret[i][j][k][l][m] = []
                        for q in range(length):
                            p[i][j][k][l][m].append(model.continuous_var(
                                lb=0, ub=1, name=f'p_{i}_{j}_{k}_{l}_{m}_{q}'))
                            p_ret[i][j][k][l][m].append(0)

    # 对偶约束
    lambda_1 = {}
    mu_1 = {}
    for i in range(length):
        lambda_1[i] = {}
        mu_1[i] = {}
        for j in range(length):
            lambda_1[i][j] = []
            mu_1[i][j] = []
            for k in range(length):
                lambda_1[i][j].append(model.continuous_var(
                    lb=0, name=f'lambda_{i}_{j}_{k}'))
                mu_1[i][j].append(model.continuous_var(
                    lb=0, name=f'mu_{i}_{j}_{k}'))
    alpha = model.continuous_var(name='Alpha')

    # 约束1️⃣
    for i in n.nodeList:
        nodeID = i.nodeID                      # --> i
        for adjacentNodeID in i.adjacentNodes:  # --> j
            dualFlowFraction = model.linear_expr()
            dualFlowFraction -= n.getLinkCapacity(nodeID, adjacentNodeID) * alpha
            for sd_p in sd:
                sd_start, sd_end = sd_p        # --> m, n
                if sd_start == sd_end:
                    continue
                dualCapacity = model.linear_expr()
                # m --> sd_start
                # n --> sd_end
                # i --> nodeID
                # j --> adjacentNodeID
                dualCapacity += lambda_1[nodeID][adjacentNodeID][sd_start] + \
                    mu_1[nodeID][adjacentNodeID][sd_end]
                for tmp_m in n.nodeList:
                    tmp_m_id = tmp_m.nodeID                 # --> u
                    for tmp_n_id in tmp_m.adjacentNodes:    # --> v
                        dualCapacity -= p[sd_start][sd_end][tmp_m_id][tmp_n_id][nodeID][adjacentNodeID]
                for t in n.getNode(sd_start).adjacentNodes:
                    dualFlowFraction += lambda_1[nodeID][adjacentNodeID][sd_start] * \
                        n.getLinkCapacity(sd_start, t)
                for s in n.getNode(sd_end).adjacentNodes:
                    dualFlowFraction += mu_1[nodeID][adjacentNodeID][sd_end] * \
                        n.getLinkCapacity(s, sd_end)
                model.add_constraint(dualCapacity >= 0)
            model.add_constraint(dualFlowFraction <= 0)

    # 约束2️⃣
    # 流约束
    # x_{ij}^{mn} 从业务量起点出发，到终点结束，沿着某一条路径，出点为1，汇点为1
    for sd_p in sd:
        sd_start, sd_end = sd_p
        if sd_start == sd_end:   # --> m,n
            continue
        for i in n.nodeList:
            nodeID = i.nodeID    # --> i
            outFlow = model.linear_expr()
            if nodeID == sd_start:
                outFlow -= 1
            if nodeID == sd_end:
                outFlow += 1
            for adjacentNodeID in i.adjacentNodes:
                if adjacentNodeID == nodeID:
                    continue
                outFlow += x[sd_start][sd_end][nodeID][adjacentNodeID]
                outFlow -= x[sd_start][sd_end][adjacentNodeID][nodeID]
            model.add_constraint(outFlow == 0)

    # 约束3️⃣
    # 修复路径约束
    # 修复路径的流，从边起点到边终点
    for sd_p in sd:
        sd_start, sd_end = sd_p                             # --> m, n
        if sd_start == sd_end:
            continue
        for iNode in n.nodeList:
            nodeI = iNode.nodeID                            # --> i
            for i in n.nodeList:
                nodeID = i.nodeID                           # --> v
                for adjacentNodeID in i.adjacentNodes:      # --> t
                    flow = model.linear_expr()
                    for tmp_out in n.getNode(nodeI).adjacentNodes:
                        flow += p[sd_start][sd_end][nodeID][adjacentNodeID][nodeI][tmp_out]
                    for tmp_in in n.getNode(nodeI).adjacentNodes:
                        flow -= p[sd_start][sd_end][nodeID][adjacentNodeID][tmp_in][nodeI]
                    if nodeI == nodeID:
                        flow -= x[sd_start][sd_end][nodeID][adjacentNodeID]
                    if nodeI == adjacentNodeID:
                        flow += x[sd_start][sd_end][nodeID][adjacentNodeID]
                    model.add_constraint(flow == 0)

    model.minimize(alpha)
    # model.export_as_lp(
    #     f"lp/model_{verticesNum}_{params.get('alpha')}_{params.get('beta')}_{params.get('seed')}.lp")

    sol = model.solve()
    if sol:
        for i in range(length):
            for j in range(length):
                for k in range(length):
                    for l in range(length):
                        if sol[x[i][j][k][l]]:
                            print(
                                f"x[{i}][{j}][{k}][{l}]: {sol[x[i][j][k][l]]}")
                            x_ret[i][j][k][l] = 1
        for i in range(length):
            for j in range(length):
                for k in range(length):
                    for l in range(length):
                        for m in range(length):
                            for q in range(length):
                                value = sol[p[i][j][k][l][m][q]]
                                # 根据具体需求设置阈值
                                if value != 0:
                                    print(
                                        f"p[{i}][{j}][{k}][{l}][{m}][{q}]: {value}")
                                    p_ret[i][j][k][l][m][q] = value
    else:
        print("求解失败")
        print(model.get_solve_status())
        print(model.get_solve_details())
        return {}, {}, 0
    return x_ret, p_ret, sol.get_objective_value()
